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Abstract 

The  authors  develop  a  model  of  steady  state  thermal  transfer  in  polymer  electrolyte  membrane  fuel  cell  stacks.  The  model  is  appropriate 
for  straight  coolant  channel  unit  cell  designs  and  considers  quantities  averaged  over  the  cross-channel  direction,  ignoring  the  impact  of  the 
gas  and  coolant  channel  geometries.  At  steady  state  and  under  the  assumption  that  the  membrane  and  electrodes  are  infinitesimally  thin,  a 
simple  description  can  be  made  of  the  temperature  distribution  in  the  cells.  The  model  provides  estimates  on  two  important  quantities:  the 
local  temperature  difference  between  coolant  and  membrane,  and  the  spread  of  heat  from  an  anomalously  hot  cell  to  its  neighbours  in  a 
stack  environment.  The  former  question  is  easy  to  address  after  the  model  is  presented.  The  latter  requires  small  argument  Laplace  transform 
asymptotics  that  correspond  to  a  large  scaled  heat  transfer  coefficient  to  the  coolant  channels.  Results  of  computational  approximation  of  the 
model  are  also  shown  and  compared  to  the  asymptotics. 

©  2005  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

A  polymer  electrolyte  membrane  fuel  cell  (PEMFC)  is  an 
electrochemical  device  in  which  the  energy  of  the  chemical 
reaction  is  converted  directly  into  electricity.  By  combining 
hydrogen  fuel  with  oxygen  from  air,  electricity  is  formed 
without  combustion  of  any  form.  Water  and  heat  are  the  only 
byproducts  when  hydrogen  is  used  as  the  fuel  source.  Further 
details  of  general  fuel  cell  operation  can  be  found  in  [1]. 

Computational  modelling  of  fuel  cell  operation  has  been 
seen  as  a  way  to  perform  design  optimization  more  efficiently 
than  by  experimental  testing  in  certain  situations.  Early  low¬ 
dimensional  averaged  models  of  unit  cell  performance  were 
developed  in  [2]  and  elsewhere.  A  modern  version  of  this  kind 
of  model  was  recently  developed  by  our  group  [3].  Three- 
dimensional  computational  models  have  been  developed  in 
[4-8],  for  example.  These  are  finite  volume  computational 
tools  that  describe  the  coupled  mass  transport  and  electro¬ 
chemistry  in  unit  cells. 


*  Corresponding  author.  Tel.:  +1  604  822  5784. 

E-mail  addresses:  kpromisl@rnath.msu.edu  (K.  Promislow), 
wetton@math.ubc.ca  (B.  Wetton). 

0378-7753/$  -  see  front  matter  ©  2005  Elsevier  B.V.  All  rights  reserved, 
doi:  10. 1016/j.jpowsour.2005.02.032 


A  few  authors  [9,10]  have  concentrated  specifically  on 
computational  models  of  heat  transfer  in  fuel  cells.  Accu¬ 
rate  approximation  of  temperature  profiles  is  of  interest  in 
the  application  since  current  fuel  cell  membranes  degrade  if 
the  temperature  increases  10-20  °C  past  standard  operating 
points.  Thermal  profiles  also  determine  (and  are  affected)  by 
where  condensation  occurs,  which  also  has  a  large  impact  on 
performance.  Despite  the  importance,  there  have  been  few  at¬ 
tempts  to  model  the  effects  of  cell-to-cell  coupling  in  a  stack 
environment.  A  first  attempt  was  made  in  our  own  work  [11]. 
The  present  paper  can  be  viewed  as  a  mathematical  analysis 
of  a  simplified  version  of  the  model  presented  there  that  gives 
insight  into  the  way  heat  is  spread  from  anomalously  hot  cells 
in  a  stack. 

The  main  simplification  of  the  model  is  that  it  considers 
heat  transfer  only  in  an  averaged  sense  over  the  cross-channel 
geometry,  as  discussed  in  Section  2.2.  This  simplification 
limits  the  model  to  qualitative  accuracy.  Another  significant 
limitation  to  the  model  is  that  it  does  not  couple  the  computed 
temperature  profiles  to  performance.  This  means  that  it  can¬ 
not  predict  phenomena  of  local  thermal  runaway,  caused  by 
such  mechanisms  as  the  feedback  loop  between  membrane 
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drying  and  heat  production.  However,  the  ability  of  the  model 
to  predict  the  spread  of  heat  in  a  stack  environment  does 
give  insight  in  to  designs  where  such  runaway  conditions 
can  occur  and  how  to  develop  designs  that  can  suppress 
them.  Currently,  there  is  no  experimental  data  to  validate 
the  model,  which  was  developed  in  part  to  allow  experi¬ 
mentalists  to  investigate  what  kind  of  measurements  would 
give  useful  insight  into  quantitative  aspects  of  heat  transfer 
in  fuel  cell  stacks.  Such  experiments  are  currently  being 
undertaken. 

In  Section  2,  we  develop  the  basic  model.  In  Section  3,  we 
nondimensionalize  the  equations  and  present  two  alternate 
reformulations  of  the  equations,  one  suitable  for  computa¬ 
tion  and  the  other  for  analysis.  In  Section  4,  we  present  the 
small  argument  Laplace  transform  analysis  that  leads  to  an 
approximate  analytic  formula  for  the  extent  of  the  spreading 
of  heat  from  an  anomalously  hot  cell  in  a  stack.  In  Section  5, 
we  compare  computational  and  analytic  solutions.  Most  of 
the  notation  and  values  of  physical  parameters  are  given  in 
Section  2. 


2.  Basic  model 

We  consider  a  unit  fuel  cell  as  shown  in  Fig.  1,  which  is 
not  drawn  to  scale.  The  electrodes  and  membrane  together 
(known  collectively  as  the  membrane  electrode  assembly  or 
MEA)  are  thinner  than  the  plates.  We  model  cells  with  straight 
(in  x)  coolant  channels  such  as  the  Ballard  Mk  9  hardware, 
although  other  designs  such  as  serpentine  channels  are  in  use 

[5]. 

Unit  cells  are  stacked  together  in  series  to  make  a  fuel  cell 
stack.  The  anode  plate  of  one  cell  is  put  against  the  cathode 
plate  of  the  next.  The  combined  plate  is  known  as  a  bipolar 
plate.  A  schematic  is  shown  in  Fig.  2.  We  will  use  subscripts 
j  to  indicate  the  cell  number  in  the  stack. 

We  make  the  following  assumptions: 
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Fig.  1.  Schematic  of  unit  fuel  cell. 
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Fig.  2.  Cross-section  (x-y)  of  stacked  unit  cells. 


1 .  The  MEAs  can  be  taken  to  be  infinitesimally  thin.  The  ef¬ 
fect  of  MEA  structure  can  be  added  as  done  in  [  1 1  ] .  These 
effects  are  not  insignificant  (the  majority  of  the  heat  is 
produced  in  the  cathode  electrode  due  to  the  ORR  over¬ 
potential).  However,  we  choose  to  neglect  these  effects  for 
simplicity  in  the  current  model  which  aims  at  qualitative 
description  of  stack  level  thermal  coupling. 

2.  Heat  is  generated  only  in  the  MEA  layer  and  that  the  heat 
flux  QJ  is  given  and  constant  in  x  and  z  for  each  unit  cell 
but  that  it  may  vary  between  cells  in  the  stack. 

3.  Heat  is  removed  only  by  the  coolant  and  that  local  (y— 
Z  plane)  conduction  to  the  coolant  is  dominant.  That  is, 
heat  conduction  in  the  x-direction  of  the  plates  can  be 
neglected.  This  will  be  true  except  in  an  asymptotically 
small  region  near  cell  inlet  and  outlet  where  corrections 
are  needed  to  match  to  given  insulating  or  radiative  con¬ 
ditions.  We  neglect  these  small  end  effects. 

4.  The  geometric  effects  of  the  gas  channels  can  be  ignored 
and  the  thermal  transport  from  electrodes  to  coolant  chan¬ 
nels  through  the  plates  can  be  described  approximately  by 
a  one-dimensional  conjugate  heat  transfer  process  with 
standard  transfer  coefficients. 

We  also  neglect  stack  end  effects  (for  example,  convection 
from  end  cell  surfaces).  These  effects  are  known  to  be  im¬ 
portant  [12],  but  our  work  here  concentrates  on  the  thermal 
coupling  effects  of  cells  away  from  stack  ends. 

Assumptions  1  and  2  above  are  made  just  for  the  con¬ 
venience  of  the  analytic  investigations  of  this  paper.  Some 
preliminary  results  of  a  more  complete  model  are  given  in 
[11],  where  the  temperature  distribution  through  the  MEA  is 
modelled  and  approximate  heat  fluxes  from  membrane  resis¬ 
tance,  condensation  and  overpotential  losses  varying  (in  x) 
are  used  from  an  existing  performance  model. 

Assumptions  3  and  4  are  what  lead  to  great  simplifications 
in  the  model  equations.  They  are  discussed  in  more  details 
in  Sections  2.1  and  2.2.  We  take  all  quantities  to  be  averaged 
over  the  cross-channel  ’-direction  and  use  T J(x)  to  denote  the 
average  coolant  temperatures,  fij{x)  for  plate  temperatures  at 
the  coolant  and  9j(x)  for  MEA  temperatures.  All  temperatures 
are  taken  relative  to  the  inlet  coolant  temperature.  We  arrive  at 
a  system  of  ordinary  differential  algebraic  equations  (DAE) 
in  these  variables  below. 
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We  use  the  following  physical  constants: 

c:  Coolant  heat  capacity  (J  kg-1  K-1). 

Qc :  Coolant  flux  per  unit  z  in  each  cell  (kg  s-1  m-1). 

A:  Coolant  transfer  factor  5600  W  m-2  K- 1  derived  in  Sec¬ 
tion  2.2. 

k\  Bipolar  plate  thermal  conductivity.  Considerable  vari¬ 
ability  in  this  value  can  be  found  for  different  plate  con¬ 
struction.  Varying  values  for  through  plane  y  and  in¬ 
plane  x-z  are  also  reported.  We  take  the  single  value 
1.29  W  m-1  K-1  [13]  that  corresponds  to  graphite  20^-0 
mesh  as  a  base  value  for  this  study. 

L\  Single  plate  thickness  (1  mm). 

Lc:  Cell  length  (1  m). 

Q :  Heat  flux  7000  Wm- 2  as  a  base  value.  This  roughly 
corresponds  to  a  current  density  of  10,000 Am-2  (or 
1  A  cm-2)  at  about  50%  overall  efficiency. 

7a:  Target  coolant  temperature  increase  from  inlet  to  outlet 
for  an  average  cell.  We  take  7a  =  15  K.  This  determines 
cQc  as  shown  below. 

The  coolant  temperature  increases  according  to  how  much 
heat  it  absorbs  locally: 


integrating  (1)  from  inlet  to  outlet  and  specifying  the  target 
temperature  increase,  we  obtain  the  value 


(5) 


for  the  combination  cQc,  which  can  be  considered  to  be  the 
thermal  carrying  capacity  of  the  coolant  flow,  per  unit  width. 
T(x)  is  linear  from  0  to  7a-  At  every  x,  the  plate  temperature 
by  the  coolant  j3  is  larger  than  the  coolant  by 


Q 

—  &  1.25  K 
A 


and  the  MEA  temperature  6  is  larger  than  this  plate  temper¬ 
ature  by 


LQ 

2k 


2.71  K 


using  the  base  parameters  above.  The  main  objective  of  this 
work  is  to  determine  how  anomalously  large  heating  in  one 
cell  in  a  stack  of  cells  at  base  conditions  will  be  distributed 
through  the  stack. 


cQc^  =  MPj-TjX  (l) 

where  the  RHS  above  describes  the  local  heat  flux  in  to  the 
coolant,  expressed  by  the  flux  from  the  edge  of  coolant  chan¬ 
nel  to  its  average  temperature.  We  now  consider  two  heat  flux 
balances,  the  first  at  the  MEA  and  the  second  at  the  coolant: 

Qj  =  fa-Pj)+lVj-P  7-1)  (2) 


2.1.  Heat  removal 

In  this  section,  we  show  that  the  heat  removal  by  ther¬ 
mal  conduction  through  the  graphite  plates  in-plane  x-z  is 
negligible  as  well  as  the  heat  removed  by  the  reactant  gas 
streams. 

We  consider  the  oxidant  stream  with  pure  air.  At  inlet,  the 
molar  flow  rate  of  oxygen  per  unit  z  is 


MPj  -  Tj)  -  5(0j  -  Pj)  +  5(dj+l  -  Pj)  (3)  Sjf, 


Here,  we  have  used  simple  expressions  for  one-dimensional 
steady  heat  conduction.  The  two  terms  on  the  right-hand  side 
of  (2)  represent  heat  flux  leaving  the  MEA  up  and  down, 
respectively.  In  (3),  the  two  terms  on  the  right-hand  side  rep¬ 
resent  the  heat  flux  flowing  from  cell  j  and  cell  j  +  1,  respec¬ 
tively,  into  coolant  channel  j. 

We  have  boundary  conditions 

Tj(  0)  =  0  (4) 

for  each  j  (recall  temperatures  are  all  relative  to  the  coolant 
inlet  temperature  common  to  all  cells  since  the  flow  is  from 
a  common  source  distributed  by  a  header).  Eqs.  ( 1 )— ( 3 )  are  a 
first-order  DAE  to  be  solved  from  inlet  x  =  0  with  conditions 
(4)  to  outlet  x  =  Lc.  In  the  theoretical  analysis  that  follows, 
we  assume  a  stack  with  an  infinite  number  of  cells.  In  com¬ 
putations,  we  can  consider  a  finite  stack,  periodic  in  j.  These 
are  conveniences  in  this  paper  where  we  are  not  explicitly 
considering  stack  end  effects. 

If  we  consider  a  uniform  stack  ( Q  /  =  Q  for  all  j),  then 
Eqs.  (l)-(3)  have  a  simple  solution  common  to  all  cells/  By 


where  S  is  the  stoichiometric  flow  rate  (the  dimensionless 
flow  rate  to  the  minimum  to  generate  the  desired  current). 
We  take  S  =  2.  I  is  the  current  density  of  the  cell  which 
we  take  to  be  10,000  Am-2  and  T  is  Faraday’s  constant 
(96485  C  mol-1).  The  heat  removed  by  the  gas  stream  at  con¬ 
stant  pressure  is  approximately 


1  J Lc 
— S— 
X0  ^ 


cn7a. 


where  XQ  ~  0.2  is  the  molar  fraction  of  oxygen  in  air 
and  cn  is  the  specific  heat  capacity  of  nitrogen  cn  = 
29.1  J  mol-1  K  1  [13].  The  heat  removed  by  the  coolant  per 
unit  width  is  QLC.  The  ratio  of  heat  removal  by  the  cathode 
gas  stream  to  the  coolant  is  then  approximately 


S  7cn7a 
~X0  QAT 


0.016. 


This  is  small  enough  to  be  neglected. 
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We  now  approximate  the  heat  transferred  per  unit  width 
by  plate  conduction,  which  can  be  found  to  be 

2Lk 

-tTa- 

The  ratio  of  this  heat  transfer  to  that  of  the  coolant  is 
2  Lk  7 

~qj2.Ta  ^  3.7  x  10-7,  (6) 

which  can  be  neglected  even  if  the  in-plane  conductivity  of 
the  plates  is  two  orders  of  magnitude  larger  than  through 
plane  which  has  been  reported  [12].  This  justifies  the  neglect 
of  the  in-plane  thermal  transport  in  the  plates  away  from  cell 
ends. 

2.2.  Averaging  over  cross-plane 


nels).  We  assume  the  flux  into  the  coolant  is  uniform 
around  its  circumference  and  that  the  coolant  flow  is  lam¬ 
inar.  Using  standard  conjugate  heat  transfer  theory  [14],  we 
find 


Q  —  =  NU—(P,  -  Ti ), 
xD  D  J  J 


where  Na  is  a  dimensionless  number  that  depends  on  the 
coolant  channel  shape,  48/1 1  for  the  case  of  the  circular  chan¬ 
nels  considered  here,  Lw  is  the  spacing  between  the  channels 
which  we  take  to  be  1  mm  and  kq  is  the  conductivity  of  the 
coolant  taken  to  be  0.41  W m_l  K-1  [12].  This  leads  to  the 
value 

A  =  Nu—  «  5600  Wm-2 K-1 


The  cross-plane  (y-z)  geometry  is  shown  in  Fig.  3.  We 
discuss  briefly  here  the  reasoning  that  leads  to  the  one¬ 
dimensional  approximation  of  the  cross-plane  transport  used 
in  the  model  equations  above.  If  the  effect  of  the  gas  channel 
geometries  are  neglected  and  the  temperatures  at  the  MEA  (0) 
and  the  temperatures  at  the  coolant  channel  edge  and  centre 
line  of  the  bipolar  plates  (ft)  can  be  taken  to  be  constant  in  z, 
we  are  led  to  formulas  (2)  and  (3).  These  assumptions  clearly 
are  not  rigorously  valid  and  are  made  for  the  convenience  of 
the  simplified  analysis  to  follow.  The  results  we  obtain  are 
limited  to  qualitative  accuracy  only. 

If  we  consider  a  net  average  thermal  flux  density  Q  coming 
from  the  MEAs  above  and  below  the  coolant  channel  shown 
in  Fig.  3,  we  see  that  the  average  thermal  flux  density  into  the 
coolant  is 


used  in  the  model  above. 

3.  Nondimensionalization  and  formulations 

We  solve  (2)  for  Of. 

=  +  +  TO 

and  use  this  result  in  (3)  to  obtain 

(  2  LA\  2  LA 

Pj+ 1  ~  (^2  +  ~ “ J  Pj  +  Pi~ i  +  yyTi 
=  ~(Qj+Qj+i).  (8) 


where  Lw  is  the  MEA  width  associated  with  one  chan¬ 
nel  (that  is,  the  cell  width  divided  by  the  number  of  chan¬ 


This  equation  with  ( 1 )  is  used  for  the  numerical  computations 
described  in  Section  5.  Values  of  the  MEA  temperatures  0] 
are  post-processed  using  (7). 

We  further  combine  (1)  with  (3)  to  eliminate  the  coolant 
temperatures  T /  and  obtain  an  equation  for  the  plate  temper¬ 
atures  /3  alone: 


Lw 


Q/2  Cathode  Electrode 


Oxidant 

Channel 

Coolant  Channel 

O  ,  Average  Temperature  T . 

Fuel 

Channel 

t 

- ^  z 

|  Q/2  Anode  Electrode 

-^{/J;+1-2(l  +  W)$j  +  $j-i} 

w 

=  KL(Pj+\  -  2 Pj  +  Pj- 1)  +  (Qj  +  Qj+ 1),  (9) 

where  the  dot  denotes  differentiation  with  respect  to  x,  and 
W  —  LA/K,a  dimensionless  ratio  of  the  temperature  increase 
through  the  plates  to  that  of  the  plates  to  the  coolant. 

We  now  scale  the  equations  and  variables.  Dimensionless 
variables  are  denoted  by  hats. 

Flux :  We  scale  the  fluxes  to  make  them  order  unity: 

Qj  =  QQj’ 


Fig.  3.  Cross-plane  (y-z)  view  of  a  unit  cell.  The  channel  structure  repeats 

in  the  ^-direction.  Unit  cells  are  stacked  in  the  ^-direction.  where  Q  is  the  representative  value  of  the  previous  section. 
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Temperature-.  The  scaling  of  temperature  is  dictated  by  the 
thermal  scaling  above 


This  can  be  seen  by  inserting  the  form  above  into  Eq.  (13)  for 
coolant  channels  not  adjacent  to  the  anomaly  (j  not  equal  to 
0  or  —1).  For  these  cells,  the  right-hand  side  of  (13)  is  zero 
and 


Length:  The  natural  length  scaling  follows: 


cQcL  „ 

x  =  - X, 

Wk 


where  the  length  scale  above  can  be  further  simplified  to 


Gj+]  -  2  ^1  +  ^  j  GJ  +  G'~ 1 1  =  0 

for  j  >  0  and  a  similar  expression  for  j  <  —  1 .  This  rela¬ 
tionship  can  be  satisfied  for  all  j  >  0  (and  the  equivalent 
relationship  for  j  <  —1)  if  G(s)  solves 


C(s)G 


1/2 


LcQ 

TaA 

with  our  parameters.  In  Section  4,  an  asymptotic  theory  for 
large  x  is  presented.  We  have  x  e  [0,  ATa/Q]  ta  [0,  12] 
for  the  base  parameters  and  numerical  experiments  show 
the  asymptotics  are  valid  at  outlet.  The  asymptotic  regime 
corresponds  to  large  ATa/Q  which  can  be  considered  to 
be  a  scaled  coolant  channel  thermal  transfer  coefficient. 

In  scaled  variables,  dropping  the  hats,  (9)  becomes 

{Pj+\  —  2(1  +  W)foj  +  Pj~ i}  +  Pj+ 1  —  2 Pj  +  Pj- 1 
=  -(Qj  +  Qj- 1).  (10) 

This  equation  is  considered  analytically  in  the  next  section. 

4.  Analysis 


9  (  Ws  \ 

G2  —  2  (l  +  - — -JG+ 1=0  (14) 

with  |G|  <  1.  Consider  (13)  at  j  =  0  (or  equivalently  at  j  = 
—  1)  to  determine  C(s): 

cWG-^(^-  2(i  +  2^-)g  +  0)  =  4 

In  the  bracket  above  add  and  subtract  1,  then  use  (14)  to 
simplify  the  expression  to 

C(s)G“1/2(G  -  1)  =  --, 


CW_  J  G-1/2_G1/2’ 


Consider  now  (8)  at  x  —  0,  where  7j(0)  =  0  for  all  j.  Let 
b  j  denote  Pj(O').  We  obtain  under  our  scaling: 

bj+i  ~  2(1  +  W)bj  +  bj—i  =  -( Qj  +  Qj+X).  (11) 

We  take  the  Laplace  transform  of  (10): 

(s  +  1)3/+1  —  2(5  +  sW  +  1  )Pj  +  (5  +  l)Pj—l 
—  bj+ 1  —  2(1  +  W)b  j  +  bj-i  —  -(Qj  +  Qj+ 1) 

=  -^+l)(Gj+Gy+t),  (12) 

where  we  have  used  tildes  for  the  Laplace  transform  and 
where  (11)  was  used.  After  dividing  by  s  +  1,  we  obtain: 

(  sW  \  -  1 

Pj+ 1  -  2  ^1  +  Pj  +  Pj- 1  =  --(Qj+  Qj+ 1). 

(13) 


and  so  the  following  is  an  explicit  form  of  the  Laplace  trans¬ 
form  for  the  plate  temperatures: 

«I)=;c^ircw^G(s)l",/2'  <15) 

For  fixed  j,  the  expression  (15)  is  analytic  except  on 
the  negative  real  axis.  It  has  branch  points  at  s  —  0  and 
s  =  —1/(2  +  W)  and  has  algebraic  growth  as  |s|  — >  oo.  It 
satisfies  the  conditions  of  Watson’s  lemma  [15]  and  the  large 
x  behaviour  is  determined  by  the  small  s  asymptotics  of  the 
transform.  We  consider  the  form  of  G  for  s  small  in  (14): 

G2  -2(1  +  Ws)G+  1  «  0, 

SO 

G^l  +  Ws  -  V(1  +  Ws)2  -  1  «  1  -  V2Ws. 


We  now  consider  an  anomalous  cell  (j  =  0)  in  an  infi¬ 
nite  stack  generating  uniformly  more  heat  than  any  other,  so 
Qo  =  1  but  Q  j  —  0  for  all  other  j.  Resulting  temperatures 
are  increases  from  base  temperatures.  Eq.  (13)  has  a  solution 
of  the  form: 

Pj(s)  =  C(s)Glj+l/21. 


Now,  (15)  becomes 


|.  ^  1  -\(j-l/2)\Q2Ws 

J  s3/2V2W 


(16) 


where  the  alternate  form  G  &  e  ^2Ws  has  been  used  in  the 
approximation  of  G^-1'2!  in  order  to  obtain  (16).  The  ex- 
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pression  above  has  a  known  inverse  transform: 


x  (ij'  -  WI) }  ■  <i7> 


where  erfc  is  the  complementary  error  function.  Since  (16) 
is  accurate  to  C)(~Js)  the  expression  (17)  has  errors  0(  1  /  -Jx) 
for  fixed  j. 

The  expression  (17)  is  analogous  to  the  solution  of  the  heat 
equation  in  a  semi-infinite  rod  forced  by  a  fixed  heat  flux  in  at 
the  end  (see  [16],  for  example),  where  x  plays  the  role  of  time 
and  j  of  a  discretized  spatial  variable.  It  is  clear  the  physical 
models  are  also  similar. 

We  note  that  in  this  setting  with  one  anomalously  heated 
cell,  plate  temperature  increases  grow  like  */x.  In  a  uniformly 
heated  stack,  base  temperatures  increase  linearly.  The  expres¬ 
sion  above  shows  that  the  anomalous  heat  spreads  out  to  a 
characteristic  number  of  adjacent  cells  J  characterized  by 


If  we  consider  the  spread  at  the  outlet  of  a  cell  of  length  Lc 
then  we  obtain,  recalling  the  length  scaling: 


J  « 


(18) 


with  our  parameters.  Note  that  this  expression  does  not  de¬ 
pend  on  A  and  increases  with  k/L.  These  predictions  are 
verified  in  the  computations  shown  in  the  next  section. 

We  remark  that  we  could  have  used  the  expression 


1  17-1/21 

s3/2vTw  s 


instead  of  (16)  leading  to  the  expression 


-17-1/21 


which  is  asymptotically  as  accurate  as  (17).  However,  this 
form  gives  negative  (unphysical)  values  for  small  x  (out  of 
the  asymptotic  range). 


5.  Computations 

We  consider  a  2M-cell  periodic  stack  using  M  computa¬ 
tional  cells  labelled  1  (anomaly)  up  to  M  using  symmetry 
conditions.  The  computations  are  based  on  an  implicit  dis¬ 
cretization  of  (1)  and  (8),  considered  as  a  linear  system  of 
2 M  unknowns  Tj  and  /i;-  at  the  next  space  step  in  x.  MEA 


Fig.  4.  Temperatures  from  an  anomalously  hot  centre  cell,  base  conditions. 
The  solid  line  is  the  anomalous  cell  j  =  1 ,  dashed  line  is  j  =  2  and  dotted 
line  is  j  =  3.  Symbols  represent  the  asymptotic  solutions. 


temperatures  6/  can  be  reconstructed  using  (7).  We  use  a  grid 
of  N  =  100  in  the  channel  length  discretization. 

As  an  example,  we  take  a  17-cell  stack  (M  =  8)  and  take 
Qq  =  Q  but  all  other  Q  j  =0.  This  simulates  a  stack  with  a 
centre  cell  running  very  hot.  On  its  own  (isolated  with  one 
coolant  channel),  it  would  run  7a  degrees  hotter  at  outlet  than 
normal.  The  temperatures  computed  from  the  study  should 
be  thought  of  as  increases  from  normal  temperature  profiles. 
In  Fig.  4,  we  show  the  results  for  the  base  parameters;  in 
Fig.  5,  the  results  if  A  is  doubled  from  the  base  case  (coolant 
thermal  conductivity  doubled);  and  in  Fig.  6,  the  results  if  k/L 
is  doubled  (plate  conductivity  doubled  or  thickness  halved). 
Note  that  doubling  A  has  little  effect  on  the  solution  spread  but 
doubling  k/L  spreads  out  the  heat  from  the  anomalous  cell  to 
more  neighbouring  cells.  This  agrees  with  the  approximate 
analytic  predictions  above. 


Fig.  5.  Temperatures  from  an  anomalously  hot  centre  cell,  doubled  A.  The 
solid  line  is  the  anomalous  cell  j  =  1 ,  dashed  line  is  j  =  2  and  dotted  line 
is  y  =  3.  Symbols  represent  the  asymptotic  solutions. 
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Fig.  6.  Temperatures  from  an  anomalously  hot  centre  cell,  doubled  k/L.  The 
solid  line  is  the  anomalous  cell  j  =  1 ,  dashed  line  is  j  =  2  and  dotted  line 
is  7  =  3.  Symbols  represent  the  asymptotic  solutions. 

In  Figs.  4-6,  the  approximate  analytic  solution  for  /I  is 
shown  for  the  equivalent  parameters.  The  agreement  is  very 
good  for  larger  x  values,  where  the  analytic  results  are  asymp¬ 
totically  valid. 

6.  Summary 

We  have  developed  models  describing  the  temperature 
distribution  in  a  fuel  cell  stack  under  several  assumptions. 
Approximate  analysis  of  the  model  leads  to  an  analytic  ex¬ 
pression  (18)  for  the  characteristic  number  of  cells  anomalous 
heat  is  spread  to. 
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